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Abstract 

In this paper, we study the two dimensional lattice Boltzmann BGK model (LBGK) 
by analytically solving a simple flow in a 2 -D channel. The flow is driven by the move- 
ment of upper boundary with vertical injection fluid at the porous boundaries. The 
velocity profile is shown to satisfy a second-order finite-difference form of the simplified 
incompressible Navier-Stokes equation. With the analysis, different boundary condi- 
tions can be studied theoretically. A momentum exchange principle is also revealed 
at the boundaries. A general boundary condition for any given velocity boundary is 
proposed based on the analysis. 
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1 Introduction 



The lattice Boltzmann equation (LBE) method has become a promising tool for simulation 
of transport phenomena in recent years. Great success has been achieved in applying LBE 



flow through porous media 
H 0, m, [0| II, reaction- 



for various flow problems such as hydrodynamics |T|, 0, |3], 
[]|, [f| 0], magnetohydro dynamics || ||, |l(|, multiphase flow 
diffusion equation (jDJ, and particle suspensions [16|. Compared to its precursor, the lattice 
gas automata (LGA), LBE method is more computationally efficient using current parallel 
computers, and some artifacts like non-Galilean invariance in LGA can be eliminated in 
LBE. Careful qualitative comparisons of LBE with traditional computational fluid dynamics 
methods showed that the method is accurate, and capable of simulating complex phenomena 
§1,0,0,0. 

Although it has been proved | 



21, 22, 18] that the LBE recovers the Navier-Stokes 



equation with a second-order of accuracy in space and time in the interior of flow domain, the 
real hydrodynamic boundary conditions have not been well studied. The no-slip condition 
on the wall is one of the examples. The classical model for the non-moving wall boundary 
condition in LBE is the so-called bounce-back rule borrowed from the LGA method. Under 
the bounce-back rule, all the particles colliding with walls bounce back to the flow domain 
in the same direction, Theoretical discussion and computational experience indicate that 
the bounce-back rule actually gives a zero velocity half way between the bounce back row 
and the first row in the flow and it introduces an error of first-order in lattice spacing for 
LGA and LBE p3|, 24, |l8|. To solve this problem, various boundary conditions fpl 



have been proposed to replace the bounce-back rule and progress has been made. In [25], 
a boundary condition for the triangular (FHP) LBGK model was proposed for any given 
velocity boundary, the boundary condition generated results of machine accuracy for plane 
Poiseuille flow. In [26], a non-slip boundary condition and prescribed pressure or velocity 
inlet condition for the 3-D 15-velocity direction LBGK model were proposed and results of 
good accuracy for various flows are achieved. However, due to the lack of the fundamental 
physical reasoning and the lack of mathematical analysis, these schemes have not revealed 
the general nature of boundary conditions. 

Recently we have developed [27] a new technique to analytically solve the LBGK equa- 
tion for 2-D Poiseuille flow and Couette flow. This technique provides us a useful tool 
to analyze the error generated by various boundary conditions. In this study, we will ex- 
tend our effort to include a flow with velocities in both directions. In addition, a general 
boundary condition for any given velocity straight boundary is proposed. 



2 Governing Equation 

In this study, we will only use the square lattice LBGK model. The procedure can be easily 
extended to the triangular lattice model, for which study is easier due to the smaller number 
of velocity directions than that of the square model. 
The model is expressed as: 

f l (x + 5e l ,t + 5)-f i (x,t) = -^[f i ( x ,t)-ft q) (x,t)}, i = 0,l,...,8, (1) 
where the equation is written in physical units. Both the time step and the lattice spacing 



2 



have the value of 5 in physical units. fi(x,t) is the density distribution function along the 
direction ej at (x, t). The particle speed ej are given by = (cos(7r(i — l)/2), sin(7r(i — 
l)/2),i = 1,2,3,4, and e* = \/2(cos(vr(i - 4 - ±)/2), sin(vr(i - 4 - |)/2),i = 5,6,7,8. Rest 
particles of type with eo = is also allowed. The right hand side represents the collision 
term and r is the single relaxation time which controls the rate of approach to equilibrium. 
The density per node, p, and the macroscopic flow velocity, u = (u, v), are defined in terms 
of the particle distribution function by 

8 8 

Y,fo = p> J2 fa = p u - ( 2 ) 

i=0 t=l 

The equilibrium distribution functions / 4 (x, t) depend only on local density and velocity 
and they can be chosen in the following form (the model d2q9 ||21{j): 

A e i) 4 n 3 i 

t Q) = ^[l + 3(e J .u) + ^(e,.u) 2 -^u-u], € = 1,2,3,4 (3) 

X 9 3 

fr ] = — p[l + 3(e 4 -u) + -( ei -u) 2 --u-u], i = 5,6,7,8. 

Assume the flow is steady and 

9u „ dv , lS 

-=0, -=0, p= const, (4) 

then /j(x, t) is only a function of y. This happens when the flow is driven by boundaries 
moving in the x-direction with injection in the y-direction from the porous boundaries (see 
Fig. 1). From Eq. (Q) we have 



f{ = £[1 + 3^ + 3^-1.5$ 

fl = ^[l + 3^_i + 3t; 2 _ 1 -1.5n 2 _ 1 ] + (l-i)^ 



fl = ^[l-Suj + ^j-LOv]] 



f{ = ^[l-3^ +1 + 3t; 2 +1 -1.5n 2 +1 ] + (l-i)/i +1 (5) 



fl = ^[1 + 3%-i + 3^_i + 3«J_! + 3^.! + 9^-1^-1] + (1 - -) fi- 



fe 
ft 



7-T— [1 - 3u 7 _i + 3v 7 _i + 3-u 2 1 + 3v 2 _ 
36r 

^-[1 - 3u j+i - 3v j+ i + 3u 2 j+1 + 3v 2 4 



+ 


(1- 


;) 


ft' 










+ 


(1- 


;) 


ft' 








+ 


(1- 


i) 








r 


+ 


(1- 




J8 









fl = ^it 1 + 3u j'+i _ 3u i+i + 3u j+i + 3v j+i ~ 9u j+ iv j+1 ] + (1 - -) fl 
where // stands for the density distribution along the direction ej at y = j5. 
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According to Eqs. (|2|, y), the x-momentum can be expressed as 

/'», = fi-fi+fi-fi-rt+fi 



2p p p 
= Y Uj + (7 + Uj+1 * + 27 W-i v i-i ~ u 3+l v 3+i\ 

+(i-^)[/r 1 -/r 1 -/^ +1 +/i +1 ] 

+(i - + pUj+1 - (ft 1 - ft 1 - ft 1 + ft 1 ) 

-(ft 1 - ft 1 + ft 1 - ft 1 )] 

= Y Uj + 67 + Uj+1 ^ + 27 ^ Uj - lVj - 1 ~ n J'+ 1 ^'+ 1 ] 

+ g(l - \)[P u i-i + P u j+i - P u j}- (6) 
which further gives us 

Uj+lVj+l ~ Uj-lVj-! _ Uj+l + - 2Uj 

26 b 2 { ' 

where v = (2r — 1)5/6 is the kinematic viscosity of the fluid The above equation is 
exactly the second-order finite-difference form of the simplified incompressible Navier-Stokes 
equation under the assumption Eq. (||) and constant pressure: 

d(uv) d 2 u 
dy dy 2 
In th y-direction, it is easy to prove that 

Vj = const. (9) 

This result is obvious for an incompressible flow under the assumption Eq.(Q). 

3 Boundary Condition 

The derivation of Eq. (^) is for the interior of the flow. The same procedure can be used to 
derive the relationship of velocities near the wall. For example, at j = 1 (near the bottom 
of the flow region), we have 

pui = fl-fl + fl-fl-ti + fl 

= + f k, + u 2 ] + f [u v - u 2 v 2 ] + (1 - -)[/° - f & - / 7 2 + /£] 

6 k>t It t 

2P P r 1 P r 1 

= — Ul + — [UQ + U 2 + — MO^O - U 2 V 2 \ 

3 6r 2r 

+(1 - l -)[ P u + P u 2 - (/° - fi - f 7 + fi) - (a 2 - fi + fi - ft)\ 

= ^-ui + 7r-[«o + ^2] + -^-[uqVq - u 2 v 2 ] + ^(1 - -)[pu + pu 2 - put] + 
6 ot It 6 t 

(l-i)p[uo-uo], (10) 
r 
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where uq = — f$ + — fe — + f$)/ P- In the above derivation, it is assumed that 
equilibrium distributions at the boundary are calculated by using uo,vo- Eq. (M) further 
gives 

u 2 v 2 -uqVq u 2 + u -2ui r-1 

26 = V P + — [«o " ^ ( n ) 

Notice that So may not be equal to uq for a specific boundary condition. Therefore, for 
any t/1 and uq ^ uq, the velocity no longer satisfies the second-order difference equation 
Eq. (]?]) at the fluid layer closest to the wall. Nevertheless, if we carefully choose a boundary 
condition to force Uq = uo, the velocities will satisfy the same second-order difference 
equation as an approximation of the Navier-Stokes equation. 
Similarly, near the top of the region, we have an equation: 

U n V n - U n - 2 V n - 2 U n + M n _2 - 2u n _i T - 1 

2$ = V #2 + —^—[Un-U n \, (12) 

where u n = (/" — f£ + f£ — f§ — f™ + f$)/p- Again, the velocities will satisfy the same 
second-order difference equation if u n = u n . 

The above restriction to Uq (as well as u n ) has a profound physical meaning. If we 
analyze the momentum exchange between a wall and its nearest fluid layer, we have (note 
that the streaming step is applied to the /j after relaxation) 

AM = /i-/i_(/0_/0) 

= (1 - i)[/ 5 ° - fg\ + f (u + 3u v ) - (/° - f 7 ) 
t or 

= fl ~ fl " (1 " -)[/s 1 " fj] ~ f («i " 3«i«i). 



6r 



Addition of the last two equations gives 



or 



2AM = AM + (l-±)[fl-f°-(fl-ri)]--f(u 1 - U o) + ^-(u 1 v 1 + u v ) 

t or It 

= (2 - — )AM - -(2 - -)(«i - u ) + ^-(uivi + u v ) + 
tot It 

(1 )p(u -u ), (13) 

r 

AM = -pv Ul U ° + ^(uivi + uo^o) + (t - l)p(«o - «o)- (14) 
2 

where (iti — uq)/5 ~ du/dy near the wall. If uq = uo, this is the statement that the mo- 
mentum exchange calculated from the exchange of /j's is equal to the momentum exchange 
carried by vertical velocity plus the viscous force. This is a momentum exchange principle. 
Obviously, for any r ^= 1, if a specific boundary conditions gives uo ^ no, it gives the wrong 
momentum exchange between a wall and its nearest fluid neighbor. In other words, in order 
to guarantee the correct momentum exchange, one must choose a boundary condition in 
which uq = uo is satisfied. 

The same conclusion can be drawn for boundary condition in the vertical direction. 
Although we have shown in the last section that Vj = const in the interior of flow domain, 
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this constant does not necessarily equal to the vertical velocity at boundaries. For instance, 
at j = 1, we have 

pvi = fl-fl + fl + fl-ti-fl 

= £-[vo + v 2 ] + ^-[vl- v 2 2 ] + (1 - -)p[% + v 2 - Vl ] 
It It t 

or 

(2r - l)(u - Vl ) + (ug - v\) + 2(r - l)(v - v ) = 0, (15) 

where vo = (f® — f® + f® + f$ — f® — f$) / P- Obviously for any r / 1, the mass conservation 
is satisfied near the boundary (v% = vq) if and only if vq — Vq = 0. 

Based on the analysis given above, we propose a boundary condition for a straight 
boundary with given velocity Uq,Vq for the square lattice. We will illustrate this at the 
bottom boundary (j = 0) as an example. After streaming, /q , /{* , f® , f% , fj , /§ are known, 
and /I , /s , /g , hence /?, need to be defined. 

• Step 1: Calculation of the density p. 

The restrictions on the boundary velocity and density discussed above state that 

P = /0 + /l + /2 + /3 + /4 + /5 + /6 + / 7 ° + /8 , (16) 
PUO = A° " / 3 ° + ft ~ fi ~ /? + / 8 °, (17) 
PV = / 2 ° - /4° + /5 + /6 - /7 - /S , (18) 

Comparison of Eqs. ( |I6[ |T§| ) gives p as: 

P = 7^[/o° + /i° + / 3 ° + 2(/ 4 ° + /r° + fi)}- (19) 
i — wo 

• Step 2: Equilibrium part of the unknown distributions. 

From the given boundary velocity and the density calculated in step 1, we can calculate 
the equilibrium part of the distributions f®, f®, /g by using Eq. (Q). 

• Step 3: Non-equilibrium part of the unknown distributions. 
Substitution of ft = f- (eq) + /P', i = 1, • • • , 8 into Eqs. (0, |TJ) gives 

/5°'-/6 ' = -(/i '-/3 '-/7 ' + /8 '), (20) 

/5 0, + /6 ' = -(/2 0, -/4 '-/ 7 '-/8 '), (21) 

where ff is the non-equilibrium distribution part of . 

Furthermore, we assume the bounce-back rule is still correct for the non-equilibrium 
part of the particle distribution normal to the boundary (in this case, f® = f± )■ 
Under this condition, the other two undefined non-equilibrium distributions can be 
uniquely determined 

/ 5 ' = / 7 ' - 0.5(/f - / 3 '), (22) 
/ 6 °' = / 8 ' + 0.5(/r-/ 3 '). (23) 



6 



The introduction of non-equilibrium part is only for the purpose of discussion of the 
method, the non-equilibrium part need not be calculated explicitly. In summary, the final 
form of the unknown distributions from steps 2,3 can be determined as 

fS = /4° + ^, (24) 
fi = / 7 °-0.5(/ 1 -/3°) + ^ + ^, (25) 
ft = /8° + 0.5(/ 1 -/3 )-^ + ^. (26) 

Once all the after-streaming distributions at the boundary are determined, the collision step 
can be easily applied to all nodes including the boundary nodes. 

The above procedures have specified a new boundary condition for any given velocity 
straight boundary. This boundary condition produces the velocity profile of the exact 



solution for the plane Poiseuille flow with forcing [27]. For the special flow with Eq. (Q), 
this boundary condition yields the correct velocity profile as the solution of the difference 
equation Eqs. (0, [ll], 12) under the condition uq = uq, and u n = u n (the formula of the 



velocity profile is given in the following section). 

4 Velocity Profile 

The governing equation for the velocity profile Eq. (^) can be solved under different tangent 
velocities at upper and bottom boundaries. For simplicity, we assume Vj = vq = const, and 
uo and u n are given. Eq. (|^) can be written as 

(2 - R)u j+ i - 4 Uj + (2 + R)uj-i = 0, (27) 

where R = vq8/u. Assuming a solution of the form 

uj = \ j , (28) 

we have a quadratic equation for A: 

(2 - R)X 2 - 4A + (2 + R) = 0, (29) 



which has two solutions: 



Ao = l, Ai = |±|. (30) 



The general solution of Eq. (pTh is: 



Uj = a\\ + 6, j = 1, • • • , n — 1, (31) 

where o, b are some constants. The corresponding difference equations near bottom and top 
boundaries, Eqs. (|ll|, p^), are used to determine a, b. These two equations can be written 
in the following way: 

(2 - R)u 2 - 4uj + (2 + R)u + e = 0, where e = U } T ^ (u - u ), (32) 

At — 1 
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(2 - R)u n - 4n n _i + (2 + i?)u n _ 2 + e n = 0, where e n = — p-(u n -u n ), (33) 

zr — 1 

Substituting the general solution in Eq. into Eqs. (|3^ , [33] ) yields a linear system of 
equations for a, 6, and solving them finally gives the solution: 

Ai — 1 Ai 1 — A"i Ai 1 — Ai £() A"i — 1 6^ . . 

^ = A^T" + "A^T n ° + "a^T^ + a^T^r- (34) 

The last two term represents the error introduced at the boundaries if uq ^ uq or u n ^ u n . 
For the boundary condition introduced in Section 3, uq = uo,u n = u n , hence eo = e n = 0. 
In the special case of uq = uq = and u n = u n = U, the solution becomes: 

Uj A, — 1 , . 

u = #31- < 35 > 

It is easy to prove that this solution is a second-order approximation of the analytical 
solution 

u e Rc y' L - 1 

where Re is the Reynolds number defined as Re= vqL/v with L being the width of the flow 
region. 



5 Discussion 

It is shown in the paper that the flow velocity from the 2D LBGK simulation in the injected 
porous boundary case satisfies a second-order difference formula as an approximation of the 
Navier-Stokes equation. This gives us a better understanding of the LBGK method. The 
velocity formula near boundaries are consistent with the velocity formula inside the flow 
domain if the boundary condition is such that the distribution functions at the boundary 
give the correct boundary velocity and thereafter the momentum exchange principle is 
satisfied. A boundary condition for any given velocity boundary for the 2D square lattice 
LBGK model can be proposed based on the analysis in this paper. This boundary condition 
can be easily extended to the 3D 15-velocity direction model. 

For the triangular lattice, the boundary condition proposed by Nobel et al. p5| gives 
the correct boundary velocity. It uses three equations: J2i=o fo = P' an d Sf=i fi e i = P u a t 
a boundary node to determine three unknowns: p and two /j's which are not defined after 
the streaming step (for example, /2,/3 at the bottom boundary). The technique does not 
have an immediate extension to the square lattice because the number of unknowns on the 
boundary in the square lattice is larger than the number of restricting equations. For the 



3D 15-velocity direction lattice, Maier et al. [26] proposed a boundary condition for solid 



boundary (can be easily extended to the case with boundaries with tangential velocities). 
First, the bounce-back is used to find the unknown /j's, hence the normal velocity is zero 
and the density and tangent velocity at the wall node can be calculated, then the /j's, 
which are on the directions pointing into the flow and have a non-zero projection on the 
wall, will be adjusted to correct the tangent velocity while keeping the normal velocity and 
density unchanged. In the case where the normal velocity is zero, the boundary condition 
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proposed in this paper is reduced to that proposed by Maier et al. In the case where the 
normal velocity is non-zero, the boundary condition in this paper can be extended to the 3D 
15-velocity direction lattice. In the present discussion, analysis of the velocity profile is also 
given and the analysis provides a framework to study any boundary condition theoretically. 
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7 Figure Caption 

Schematic plot of flow driven by the movement of upper boundary with vertical injection 
fluid at the porous boundaries. Also included is a 9-bit square lattice used in this paper. 
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